自己紹介

酸素摂取量(Kuehl, 2000)

酸素摂取量(Kuehl, 2000)

• 年齢と運動療法の組み合わせごとに集団の平均と分散を推定する？
  • 例: 22歳のランニンググループとエアロビクスグループで推定
    
    
• ✖️参加者が一人しかいない年齢もあり
  • →集団の分散に関する情報は不十分
• ✖️データが存在しない年齢と運動療法の組み合わせも無数にある

解決策: 線形回帰モデル(linear regression model)を使う

• 条件付き分布\(p(y|x)\)が\(x\)の関数として滑らかに変化すると仮定

• 得られた\(x\)のデータから他の値の情報を得る

• 条件付き平均\(E[Y|x]\)はパラメータに関して線形であると定める

\[ \int yp(y|\boldsymbol{x})dy = E[Y|\boldsymbol{x}]=\beta_1x_1+...+\beta_px_p=\boldsymbol{\beta}^T\boldsymbol{x} \]

今回のモデルの仮定

• \(p(y|\boldsymbol{x})\)は年齢と酸素摂取量は線形関係
• 運動療法のグループごとに異なる線形関係を想定(異なる関係でもいい)

\[ Y_i = \beta_1x_{i,1}+\beta_2x_{i,2}+\beta_3x_{i,3}+\beta_4x_{i,4}+\epsilon_i\tag{9.1} \]

今回のモデルの仮定

今回のモデルでの\(Y\)の条件付き期待値は、\(x_{i,2}\)のとりうる値によって次のようになる

\[ E[Y|\boldsymbol{x}] = \beta_1+\beta_3\times年齢(x_2=0の場合)\\ E[Y|\boldsymbol{x}] = (\beta_1+\beta_2)+(\beta_3+\beta_4)\times年齢(x_2=1の場合) \]

年齢との線形関係は運動療法のグループ間で切片と傾きの違いがあることが仮定

酸素摂取量に対する四つのモデルの回帰直線

正規線形回帰モデル(normal linear regression model)

• \(E[Y|\boldsymbol{x}]\)が線形である
• 平均まわりでの標本のばらつきが
• 独立かつ同一の正規分布に従う(i.i.d., p.30参照)

\[ \epsilon_1,...,\epsilon_n\sim \mathrm{i.i.d. normal}(0,\sigma^2)\\ Y_i=\boldsymbol{\beta}^T\boldsymbol{x}_i+\epsilon_i \]

→\(\boldsymbol{x}_i,\boldsymbol{\beta},\sigma^2\)で条件づけたもとで観測データ\(y_1,...y_n\)の同時分布を完全に特定する

正規線形回帰モデル(normal linear regression model)

同時確率密度は式(9.2)で書ける

\[ p(y_1,...,y_n|\boldsymbol{x}_1,...,\boldsymbol{x}_n,\boldsymbol{\beta},\sigma^2)\tag{9.2}\\ \]

\[ =\Pi_{i=1}^n p(y_i|\boldsymbol{x}_i,\boldsymbol{\beta},\sigma^2)\\ =(2\pi\sigma^2)^{-n/2}\mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2\}\tag{9.3} \]

正規線形回帰モデル(normal linear regression model)

この同時確率密度は多変量正規分布を用いて書ける

\[ \{\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2\}\sim \mathrm{multivariate\ normal}(\boldsymbol{X\beta},\sigma^2\mathrm{\boldsymbol{I}}) \]

• \(\boldsymbol{y} = (y_1,...,y_n)^T\)
• \(\boldsymbol{X}\) : \(n\times p\)行列. 第\(i\)行目が\(x_i\)
• \(\mathrm{\boldsymbol{I}}\) : \(n\times n\)単位行列

正規線形回帰モデル(normal linear regression model)

\(\boldsymbol{X\beta}\)は以下で示せる

\[ \boldsymbol{X\beta}＝\begin{pmatrix}x_1\rightarrow \\ x_2 \rightarrow\\ \vdots \\ x_n \rightarrow\end{pmatrix} \begin{pmatrix}\beta_1 \\ \beta_2\\ \vdots \\ \beta_p\end{pmatrix} =\begin{pmatrix}\beta_1x_{1,1}+\dots\beta_px_{1,p} \\ \vdots \\ \beta_px_{n,1}\dots\beta_px_{n,p}\end{pmatrix} =\begin{pmatrix}\mathrm{E}[Y_1|\boldsymbol{\beta},\boldsymbol{x}_1]\\ \vdots \\ \mathrm{E}[Y_n|\boldsymbol{\beta},\boldsymbol{x}_n]\end{pmatrix} \]

\[ (2\pi\sigma^2)^{-n/2}\mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2\}\tag{9.3} \]

• 式(9.3)の密度は残差\((y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)\)を通じて\(\boldsymbol{\beta}\)に依存している

• 観測されたデータを所与とすると、残差平方和\(\mathrm{SSR(\boldsymbol{\beta})}=\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2\)を最小にすることで尤度が最大になる

• 残差平方和を最小にするには微分する

\(\hat{\boldsymbol{\beta}}_{ols}\)を求める

\(\hat{\boldsymbol{\beta}}_{ols}\) : 最小二乗推定量

• \(\hat{\boldsymbol{\beta}}_{ols}\)という値は、\(\boldsymbol{\beta}\)の「最小二乗」(ordinary least squares, OLS)推定量と呼ばれる。
  • \(\hat{\boldsymbol{\beta}}_{ols}=(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}\)
• この値は、\((\boldsymbol{X}^T\boldsymbol{X})^{-1}\)が存在するなら一意に定まる

9.1.1 酸素摂取量データに対する最小二乗推定

式(9.1)のモデルにおいて最小二乗推定量を求め、運動療法の違いを評価する

• 測定された酸素摂取量の変化: \(\boldsymbol{y}\)

\[ \boldsymbol{y}= (-0.87,-10.74,-3.27,-1.97,7.50,-7.25,17.05,4.96,10.40,11.05) \]

• 式(9.1)に定義されているベクトル\(x_1,x_2,x_3,x_4\)から行列\(\boldsymbol{X}\)を構成する

      [,1] [,2] [,3] [,4]
 [1,]    1    0   23    0
 [2,]    1    0   22    0
 [3,]    1    0   22    0
 [4,]    1    0   25    0
 [5,]    1    0   27    0
 [6,]    1    0   20    0
 [7,]    1    1   31   31
 [8,]    1    1   23   23
 [9,]    1    1   27   27
[10,]    1    1   28   28
[11,]    1    1   22   22
[12,]    1    1   24   24

9.1.1 酸素摂取量データに対する最小二乗推定

t(X)%*%X

     [,1] [,2] [,3] [,4]
[1,]   12    6  294  155
[2,]    6    6  155  155
[3,]  294  155 7314 4063
[4,]  155  155 4063 4063

y<-c(-0.87,-10.74,-3.27,-1.97,7.50,-7.25,17.05,4.96,10.40,11.05,0.26,2.51)
t(X)%*%y

        [,1]
[1,]   29.63
[2,]   46.23
[3,]  978.81
[4,] 1298.79

\(\hat{\boldsymbol{\beta}}_{ols}=(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}\)を計算する

beta.ols<- solve(t(X)%*%X)%*%t(X)%*%y
print(beta.ols)

            [,1]
[1,] -51.2939459
[2,]  13.1070904
[3,]   2.0947027
[4,]  -0.3182438

※正しい回帰直線はスライド11の(\(\beta_2\neq0,\ \beta_4\neq0\))

9.1.1 酸素摂取量データに対する最小二乗推定

\(\sigma^2\)の不偏推定量は\(\mathrm{SSR}(\hat{\boldsymbol{\beta}}_{ols})/(n-p)\)

SSR_beta.ols <- t(y)%*%y -(2*t(beta.ols)%*%t(X)%*%y) +(t(beta.ols)%*%t(X)%*%X%*%beta.ols)
n<-length(y)
p<-dim(X)[2]
print(SSR_beta.ols/(n-p))

         [,1]
[1,] 8.542477

9.2.1 準共役事前分布

\(\boldsymbol{x}_i,\boldsymbol{\beta},\sigma^2\)で条件づけたもとで観測データ\(y_1,...y_n\)の同時分布を式(9.3)に示した

\[ p(y_1,...,y_n|\boldsymbol{x}_1,...,\boldsymbol{x}_n,\boldsymbol{\beta},\sigma^2)\\ =(2\pi\sigma^2)^{-n/2}\mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2\}\tag{9.3} \]

このデータ密度は、\(\boldsymbol{\beta}\)の関数として次の様にかける

\[ p(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2)\propto \mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n\mathrm{SSR}(\boldsymbol{\beta})\}\\ =\mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n [\boldsymbol{y}^T\boldsymbol{y}-2\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{y}+\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}]\} \]

• 指数の中で\(\boldsymbol{\beta}\)は\(\boldsymbol{y}\)と同じような役割をしている
• また、\(\boldsymbol{y}\)の分布は多変量正規分布
• つまり、\(\boldsymbol{\beta}\)について多変量正規分布は共役になることが示唆

もし\(\boldsymbol{\beta}\sim \mathrm{multivariate normal}(\boldsymbol{\beta}_0,\boldsymbol{\sum}_0)\)ならば

\[ p(\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2) \propto p(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2)\times p(\boldsymbol{\beta})\\ \propto \mathrm{exp}\{-\frac{1}{2}(-2\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{y}/\sigma^2+\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}/\sigma^2) -\frac{1}{2}(-2\boldsymbol{\beta}^T\sum_0^{-1}\boldsymbol{\beta}+\boldsymbol{\beta}^T\sum_0^{-1}\boldsymbol{\beta})\}\\ =\mathrm{exp}\{ \boldsymbol{\beta}^T(\sum_0^{-1}\boldsymbol{\beta}_0+\boldsymbol{X}^T\boldsymbol{y}/\sigma^2) -\frac{1}{2}\boldsymbol{\beta}^T(\sum_0^{-1}+\boldsymbol{X}^T\boldsymbol{X})/\sigma^2)\boldsymbol{\beta}\} \]

これは、多変量正規分布の密度関数であり、平均と分散は以下で示される

\[ \mathrm{Var}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2]=(\sum_0^{-1}+\boldsymbol{X}^T\boldsymbol{X}/\sigma^2)^{-1}\tag{9.4} \]

\[ \mathrm{E}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2]=(\sum_0^{-1}+\boldsymbol{X}^T\boldsymbol{X}/\sigma^2)^{-1}(\sum_0^{-1}\boldsymbol{\beta}_0+\boldsymbol{X}^T\boldsymbol{y}/\sigma^2)\tag{9.5} \]

• 精度行列\(\sum_0^{-1}\)の要素の絶対値が小さいとき、条件付き期待値\(\mathrm{E}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2]\)は近似的に最小二乗推定量\((\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}\)に等しい
• 観測の精度(\(\sigma^2\)が非常に大きい)が非常に小さいときには、近似的に事前の期待値\(\boldsymbol{\beta}_0\)に等しくなる

\(\sigma^2\)の準共役事前分布

ほとんどの正規モデルでは、\(\sigma^2\)の準共役事前分布は逆ガンマ分布である

(5. 正規モデル参照)

\(\gamma=1/\sigma^2\)を観測の精度とすると、\(\gamma\sim\mathrm{gamma}(\nu_0/2,\nu_0\sigma_o^2/2)\)ならば

\[ p(\boldsymbol{\gamma}|\boldsymbol{y},\boldsymbol{X},\boldsymbol{\beta})\propto p(\boldsymbol{\gamma})p(\boldsymbol{y}|\boldsymbol{\gamma},\boldsymbol{X},\boldsymbol{\beta})\\ \propto[\gamma^{\nu_0/2-1}\mathrm{exp}(-\gamma\times\nu_0\sigma_0/2)] \times[\gamma^{n/2}\mathrm{exp}(-\gamma\times\mathrm{SSR}(\boldsymbol{\beta})/2)]\\ =\gamma^{(\nu_0+n)/2-1}\mathrm{exp}(-\gamma[\nu_0\sigma_o^2+\mathrm{SSR}(\boldsymbol{\beta})]/2 \]

これはガンマ分布の密度だとみなせるので、以下の結果が得られる

\[ \{\sigma^2|\boldsymbol{y},\boldsymbol{X},\boldsymbol{\beta}\}\sim \mathrm{inverse\ gamma}([\nu_0+n]/2, [\nu_0\sigma_0^2+\mathrm{SSR}(\boldsymbol{\beta})]/2) \]

ギブスサンプラーによる事後分布近似

1. \(\boldsymbol{\beta}\)の更新
   1. \(\mathrm{V}=\mathrm{Var}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^{2(s)}]\)と\(\mathrm{m}=\mathrm{E}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2{2(s)}]\)を計算する
   2. \(\boldsymbol{\beta}^{(s+1)}\sim\mathrm{multivariate normal(\mathrm{m},\mathrm{V})}\)を生成
2. \(\sigma^2\)の更新
   1. \(\mathrm{SSR}(\boldsymbol{\beta}^{(s+1)})\)を計算する
   2. \(\sigma^{2_{(s+1)}}\sim \mathrm{inverse\ gamma}([\nu_0+n]/2, [\nu_0\sigma_0^2+\mathrm{SSR}(\boldsymbol{\beta^{(s+1)}})]/2)\)を生成

9.2.2 規定事前分布と弱情報事前分布

• 回帰モデルのベイズ分析では事前分布のパラメータ(\(\boldsymbol{\beta}_0,\boldsymbol{\sum}_0\))と(\(\nu_0,\sigma^2_0\))を特定する必要がある
• これらの値を上手く見つけて、実際の事前の情報を表すのは難しい場合もある
  • 結構ありそう
• 今回の様な比較的単純なモデルでも特定する必要がある値はそこそこある
  • 相関パラメータはpが増えるに従って劇的に増加

※どんな感じで事前分布のパラメータの値を決めてるかはp.173参照

• 事前の情報が正確でない、共役事前分布で容易に表現できない
  • →他の基準によって事前分布を正当化することがある

どうする？

• 研究対象について事前の知識を全く持ってないという信念を事前分布で表現
  • 「より客観的」な結論を得る

弱情報事前分布の例として単位情報事前分布(unit information prior)をあげる

(Kass and Wasserman, 1995)

単位情報事前分布(unit information prior)

• 1個の観測地に含まれる情報と同じ量の情報を含む事前分布

  • \((\boldsymbol{X}^T\boldsymbol{X})/\sigma^2\): n個の観測値に含まれる情報の量(\(\hat{\boldsymbol{\beta}}_{ols}\)の精度)
  • \((\boldsymbol{X}^T\boldsymbol{X})/(n\sigma^2)\): 1個の観測値に含まれる情報

• よって、単位情報事前分布では\(\sum_0^{-1}=(\boldsymbol{X}^T\boldsymbol{X})/(n\sigma^2)\)

• \(\boldsymbol{\beta}_0=\hat{\boldsymbol{\beta}}_{ols}\)と定めることで、\(\boldsymbol{\beta}\)の事前分布の中心を最小二乗推定量にすることが提案されている

• 単位情報事前分布は、\(y\)に含まれる情報を少ししか使ってない

  • →バイアスがなく情報の少ない事前分布といえる

• 同様に\(\nu_o=1,\sigma_o^2=\hat{\sigma}_{ols}^2\)と定めることで、\(\sigma^2\)の事前分布の中心を\(\hat{\sigma}_{ols}^2\)とする

不変原理

• パラメータ推定は説明変数の単位の返還に関して不変であるべき

• 例: 酸素摂取量を分析するときに用いた年齢(\(x_{i,3}\))が月齢(\(\tilde{x}_{i,3}\))で事後分布は変わらないようにする

  • 年齢の事後分布と月齢の事後分布は同じはずである

• 与えられた説明変数\(\boldsymbol{X}\)について、ある\(p\times p\)行列\(\boldsymbol{H}\)があって\(\tilde{\boldsymbol{X}}=\boldsymbol{XH}\)

• 不変原理において、\(y\)と\(\boldsymbol{X}\)から\(\tilde{\boldsymbol{\beta}}\)の事後分布をえたとき、\(\boldsymbol{\beta}=\boldsymbol{H\tilde{\beta}}\)

• この条件が成立するのは、任意の正の値\(k\)に対して\(\boldsymbol{\beta}_0=0\)かつ\(\sum_0=k(\boldsymbol{X}^T\boldsymbol{X})^{-1}\)となるとき

• \(k\)の特定の仕方として、代表的なものは誤差の分散\(\sigma^2\)に関係させて\(k=g\sigma^2\)とすること

  • これはg事前分布(Zellner, 1986)の一種
  • この後のモデル選択にも用いられる

g事前分布

• g事前分布のもとで\((\boldsymbol{y},\boldsymbol{X},\sigma^2)\)を所与としたときの\(\boldsymbol{\beta}\)の条件付き分布は多変量正規分布
• 式(9.4)(9.5)は以下の様な単純な形に帰着する

\[ \mathrm{Var}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2] =[\boldsymbol{X}^T\boldsymbol{X}/(g\sigma^2)+\boldsymbol{X}^T\boldsymbol{X}/\sigma^2]^{-1}\\ =\frac{g}{g+1}\sigma^2(\boldsymbol{X}^T\boldsymbol{X})^{-1}\tag{9.6} \]

\[ \mathrm{E}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2] =[\boldsymbol{X}^T\boldsymbol{X}/(g\sigma^2)+\boldsymbol{X}^T\boldsymbol{X}/\sigma^2]^{-1}\boldsymbol{X}^T\boldsymbol{y}/\sigma^2\\ =\frac{g}{g+1}\sigma^2(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}\tag{9.7} \]

• \(p(\sigma^2|\boldsymbol{y},\boldsymbol{X})\)は逆ガンマ分布
  • \((\sigma^2,\boldsymbol{\beta})\)を事後分布から直接サンプリングすることが可能

\(p(\sigma^2|\boldsymbol{y},\boldsymbol{X})\)の導出

• \(\sigma^2\)の周辺事後密度は\(p(\sigma^2)\times(\boldsymbol{y}|\boldsymbol{X},\sigma^2)\)​に比例する
• 周辺確率の法則を用いて、積の第2項は以下の積分で書ける

\[ (\boldsymbol{y}|\boldsymbol{X},\sigma^2) =\int(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2)(\boldsymbol{\beta}|\boldsymbol{X},\sigma^2)d\boldsymbol{\beta} \]

積分の中にある二つの密度は式(9.6)を\(V\)と式(9.7)を\(m\)といて以下のように書ける

(式の展開はp.175~176参照)

\[ (\boldsymbol{y}|\boldsymbol{X},\sigma^2) =\int(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2)(\boldsymbol{\beta}|\boldsymbol{X},\sigma^2)d\boldsymbol{\beta}\\ =[(2\pi\sigma^2)^{-n/2}\mathrm{exp}(-\frac{1}{2\sigma^2}\boldsymbol{y}^T\boldsymbol{y})] \times [(1+g)^{-p/2}\mathrm{exp}(\frac{1}{2}m^TV^{-1}m)] \]

ここで、指数の項を組み合わせることで

\[ (\boldsymbol{y}|\boldsymbol{X},\sigma^2)=(2\pi)^{-2/n}(1+g)^{-p/2}(\sigma^2)^{-n/2}\mathrm{exp}(\frac{1}{2}\mathrm{SSR}_g) \]

\[ \mathrm{SSR}_g =\boldsymbol{y}^T\boldsymbol{y}-\sigma^2\boldsymbol{m}^T\boldsymbol{V}^{-1}\boldsymbol{m} =\boldsymbol{y}^T(\mathrm{I}-\frac{g}{g+1}\boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T)\boldsymbol{y} \]

• \(\mathrm{SSR}_g\)は、\(g\rightarrow\infty\)の時に減少して\(\mathrm{SSR}_{ols}=\sum(y_i-\boldsymbol{\hat{\beta}}_{ols}x_i)^2\)となる
• \(g\)の値は回帰係数の大きさを縮小し、データの過適合を防ぐ効果を持つ

• \((\boldsymbol{y}|\boldsymbol{X},\sigma^2)\)があきらかになったので、事前分布\(p(\sigma^2)\)をかけることで\(p(\sigma^2|\boldsymbol{y},\boldsymbol{X})\)は導出される
• \(\gamma=1/\sigma^2\sim\mathrm{gamma}(\nu_0/2, \nu_0\sigma_0^2/2)\)とおく
  • \(p(\sigma^2|\boldsymbol{y},\boldsymbol{X})\propto p(\sigma^2)(\boldsymbol{y}|\boldsymbol{X},\sigma^2)\)→\(p(\gamma|\boldsymbol{y},\boldsymbol{X})\propto p(\gamma)(\boldsymbol{y}|\boldsymbol{X},\gamma)\)

\[ p(\gamma|\boldsymbol{y},\boldsymbol{X})\propto p(\gamma)(\boldsymbol{y}|\boldsymbol{X},\gamma)\\ \propto \mathrm{dgamma}(\gamma,[\nu_0+n]/2,[\nu_0,\sigma_o^2+\mathrm{SSR}_g]/2) \]

\[ \{\sigma^2|\boldsymbol{y},\boldsymbol{X}\}\sim \mathrm{inverse\ gamma}([\nu_0+n]/2,[\nu_0,\sigma_o^2+\mathrm{SSR}_g]/2) \]

酸素摂取量のベイズ分析

  Xintercept     Xaerobic         Xage Xaerobic.age 
 -47.3482578   12.0988527    1.9335717   -0.2937635 

モデル選択

• 回帰分析では膨大な数の説明変数を扱うこともある
• しかし\(Y\)と本当に関連しているかは怪しい
  • 全て入れると統計分析の性能は低くなる
• \(y\)との関連が十分に確からしい変数だけをモデルに含めたい

例: 糖尿病

• 442人の糖尿病患者に関する10種類の変数\(x_1,\dots,x_{10}\)
• 参照時点から1年後にどれだけ病状が進行したのかの指標\(y\)
• 参照時点から\(y\)を予測するモデルを作成したい

モデルどうする？

• \(y\)と\(x_j\)の関係は線形ではないかもしれない
• \(x_j^2\)や\(x_jx_k\)の様な項があったほうがいいかもしれない

僕が考えた最強のモデル

• 説明変数64個

  • 主効果項: 説明変数10個

  • 交互作用項: \(\tbinom{10}{2}=45\)個

  • 二次の項: 9個(性別除く)

• 変数は標準化(\(y,\boldsymbol{X}\)が平均0,分散1)

モデル評価(方法)

• 442人の糖尿病患者をランダムに分割
  • 342人の訓練データ: \(y,\boldsymbol{X}\)
  • 100人のテストデータ: \(y_{test},\boldsymbol{X}_{test}\)
• 訓練データから推定された回帰係数を用いて、\(\hat{y}_{test}=\boldsymbol{X}_{test}\hat{\boldsymbol{\beta}}\)を生成
• \(\hat{y}_{test}\)と\(y_{test}\)を比較して性能評価

モデル評価(結果)

[1] 0.669
[1] "平均二乗誤差: 0.668828410609429"

モデル選択手順(model selection procedure)

• 回帰係数の真の値が0でないという証拠をt統計量から評価する方法がある
  • t統計量は最小二乗推定量を標準偏差でわることで得られる
  • \(t_j=\hat{\beta}_j/[\hat{\sigma}^2(\boldsymbol{X}^T\boldsymbol{X}_{j,j}^{-1})]^{1/2}\)

1. 推定値\(\hat{\boldsymbol{\beta}}_{ols}=(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}\)とそのt統計量を得る
2. もし\(j\)番目の説明変数について\(|t_j|<t_{cutoff}\)となるならば
   1. \(|t_j|\)が最小となる説明変数の番号を\(j_{min}\)として、\(\boldsymbol{X}\)から\(j_{min}\)列目を取り除く
   2. 1に戻る
3. モデルに残った全ての\(j\)についてカットオフを上回ったら終了

この方法は、変数減少法(backward elimination)の一種である

糖尿病データの場合

[1] 0.6594587

変数減少法の欠点

仮に、すべての説明変数が\(Y\)と独立な場合、この手順によってどのようなモデルが得られるのか

実際には変数間は独立だが、変数減少法によって誤った示唆がなされてしまうことがある(Berk, 1978)

ベイズ的モデル選択

• あるモデルの多くの回帰係数が0に等しいかもしれない
  • その信念を表現する事前分布を考えればよい
• 例えば、\(z\in[0,1]\)かつ\(b_j\)を正の実数として、\(j\)番目の回帰係数を\(\beta_j=z_j\times b_j\)と書く
  • \(z\)=1の時、回帰係数が0ではない

この時、回帰式は

\[ y_i=z_ib_1x_{i,1}+,\dots,+z_pb_px_{i,p}+\epsilon_i \]

酸素摂取量データでは

\[ E[Y|x,b,z=(1,0,1,0)]=b_1x_1+b_3x_3\\ E[Y|x,b,z=(1,1,0,0)]=b_1x_1+b_2x_2 \]

• \(z=(z_1,\dots,z_p)\)の値ごとに異なるモデルが対応する
• ベイズ的モデル選択では\(z\)に関する事前分布を得ることで進行する
  • \(\{\boldsymbol{z},\boldsymbol{\beta},\boldsymbol{\sigma^2}\}\)の同時事前分布

→g事前分布が使えそう

• モデルの事前分布\(p(z)\)を所与として、各回帰モデルの事後確率は次の様に計算できる

\[ p(z|y,\boldsymbol{X}) =\frac{p(z)p(y|\boldsymbol{X},z)}{\sum_{\tilde{z}}p(\tilde{z})p(y|\boldsymbol{X},\tilde{z})} \]

任意の二つのモデル比較は、事後オッズを比較することでも可能

\[ \mathrm{odds}(z_a,z_b|y,\boldsymbol{X}) =\frac{p(z_a|y,\boldsymbol{X})}{p(z_b|y,\boldsymbol{X})} =\frac{p(z_a)}{p(z_b)} \times\frac{p(y|\boldsymbol{X},z_a)}{p(y|\boldsymbol{X},z_b)} \]

• 事後オッズは、事前オッズとベイズファクターの積
• ベイズファクターは、データに対してモデル\(z_a\)が\(z_b\)に比べてどれだけ利するか解釈可能
• 事後確率\(p(y|\boldsymbol{X},z)\)を計算する必要がある

周辺確率の計算

\[ p(y|\boldsymbol{X},z)=\int\int p(y,\boldsymbol{\beta},\sigma^2|\boldsymbol{X},z)d\beta d\sigma^2\\ =\int\int p(y|\boldsymbol{\beta},\boldsymbol{X},z,\sigma^2)p(\boldsymbol{\beta}|\boldsymbol{X},z,\sigma^2)p(\sigma^2) d\beta d\sigma^2 \tag{9.8} \]

• g事前分布の一種を用いることで簡単に表現
• \(z\)を所与とし、\(p_z\)個の要素が0でないとして、\(z_j=1\)となるような変数からなる\(n\times p_z\)行列を\(\boldsymbol{X}_z\)
• \(\boldsymbol{\beta}_z\)は\(p_z\times 1\)ベクトル

• \(\boldsymbol{\beta}\)の変形g事前分布は、\(z_j=0\)となる\(j\)については\(\beta=0\)とし、かつ

\[ \{\boldsymbol{\beta}_z|\boldsymbol{X}_z,\sigma^2\} \sim\mathrm{multivariate\ normal}(\boldsymbol{0},g\sigma^2[\boldsymbol{X}_z^T\boldsymbol{X}_z]^{-1}) \]

式(9.8)で\(\boldsymbol{\beta}\)を積分

\[ p(y|\boldsymbol{X},z)=\int(\int p(y|\boldsymbol{X},z,\sigma^2,\boldsymbol{\beta})p(\boldsymbol{\beta}|\boldsymbol{X},z,\sigma^2)d\boldsymbol{\beta})p(\sigma^2)d\sigma^2\\ =\int p(y|\boldsymbol{X},z,\sigma^2)p(\sigma^2)d\sigma^2 \]

• 周辺確率\(p(y|\boldsymbol{X},z,\sigma^2)\)の計算はg事前分布でやった
• \(\gamma=1/\sigma^2\)とし、ガンマ分布の密度\(p(\gamma)\)をパラメータ\((\nu_0/2,\nu_0\sigma^2_0/2)\)を計算する

• \((\boldsymbol{X},z)\)を所与とした\((y,\gamma)\)の条件付き密度は次のようになる

\[ p(y|\boldsymbol{X},z,\gamma)\times p(\gamma)\\ =(2\pi)^{-n/2}(1+g)^{-p_z/2}\times[\gamma^{n/2}e^{-\gamma\mathrm{SSR}_g^z/2}] \times(\nu_0\sigma_0^2/2)\Gamma(\nu_0/2)^{-1}[\gamma^{\nu_0/2-1}e^{-\gamma\nu_0\sigma_0^2/2}]\tag{9.9} \]

• \(\mathrm{SSR}_g^z\)の説明変数行列は\(\boldsymbol{X}_z\)である
• 式(9.9)では\(\gamma\)に依存する箇所はガンマ分布の密度に比例

\[ \gamma^{(\nu_0+n)/2-1}\mathrm{exp}(-\gamma\times(\nu_0\sigma^2_0+\mathrm{SSR}_g^z)/2) =\frac{\Gamma([\nu_0+n]/2)}{([\nu_0\sigma^2_0+\mathrm{SSR}_g^z]/2)^{(\nu_0+n)/2-1}} \]

• 右辺にガンマ分布の密度の項もあるが、積分して消える(p.184参照)

よって\(p(y|\boldsymbol{X},z)\)は

\[ p(y|\boldsymbol{X},z)=\pi^{-n/2}\frac{\Gamma([\nu_0+n]/2)}{\Gamma(\nu_0/2)}(1+g)^{-p_z/2}\frac{(\nu_0\sigma_0^2)^{\nu_0/2}}{(\nu_0\sigma^2_0+\mathrm{SSR}_g^z)^{(\nu_0+n)/2}} \]

任意の二つのモデル\(z_a\)と\(z_b\)の比(ベイズファクター)は

\[ \frac{p(y|\boldsymbol{X},z_a)}{p(y|\boldsymbol{X},z_b)} =(1+n)^{p_{z_b}-p_{z_a}/2}(\frac{s_{z_a}^2}{s_{z_b}^2})^{1/2} \times(\frac{s_{z_b}^2+\mathrm{SSR}_g^{z_b}}{s_{z_a}^2+\mathrm{SSR}_g^{z_a}})^{(n+1)/2} \]

酸素摂取量の例

酸素摂取量のデータに対する回帰モデルは以下の通りであった

\[ E[Y_i|\boldsymbol{\beta},x_i]=\beta_1x_{i,1}+\beta_2x_{i,2}+\beta_3x_{i,3}+\beta_4x_{i,4} \]

• 運動療法の効果の有無は\(\beta_2,\beta_4\)が0かどうか
• 前節では、効果はないことが示唆

• \(\beta\)に対してg事前分布

• \(\sigma^2\)に対して単位情報事前分布

• 5つのモデルの中で最も確率が高い者は\(z=(1,1,1,0)\)

  • グループで切片が異なるが、傾きへの交互作用はないモデル

• 年齢を含む三つのモデルの事後確率を足すとほぼ1になる

  • 年齢に関する証拠は強い

• 運動療法の事後確率は0.82

  • 事前確率0.60なのでそれよりは高い確率

9.3.2 ギブスサンプリングとモデル平均

• pこの回帰係数が0かどうか考察するときを考える
  • モデルの数は\(2^p\)個
  • pが大きいと、それぞれのモデルの周辺確率を計算するのは非現実的
• 糖尿病データはp=64
  • モデルの合計は\(2^{64}\)
• ギブスサンプリングを使えば、確率が比較的高いモデルのリストや\(\beta\)のよい推定値が得られる

ギブスサンプリング

• 現在の値\(\boldsymbol{z}=(z_1,\dots,z_p)\)を所与として、新たな\(z_j\)の値を\(p(z_j|y,\boldsymbol{X},z_{-j})\)
• 完全条件付確率は\(o_j/(1+o_j)\)と書ける
• \(o_j\)は\(z_j\)が1となる条件付きオッズ

\[ o_j=\frac{\mathrm{Pr}(z_j=1|y,\boldsymbol{X},z_{-j})}{\mathrm{Pr}(z_j=0|y,\boldsymbol{X},z_{-j})} =\frac{\mathrm{Pr}(z_j=1)}{\mathrm{Pr}(z_j=1)} \times \frac{p(y|\boldsymbol{X},z_{-j},z_j=1)}{p(y|\boldsymbol{X},z_{-j},z_j=0)} \]

• \(\mathrm{\beta},\sigma^2\)も事後分布からサンプリングする必要がある
  • 9.2節の結果から直接サンプリング可能

1. \(z=z^{(s)}\)と定める
2. ランダムな順番で\(j\in\{1.,\dots,p\}\)に対し、\(z_j\)を\(p(z_j|z_{-j},y,\boldsymbol{X})\)からの標本で置き換える
3. \(z^{s+1}=z\)と定める
4. \(\sigma^{2(s+1)}\sim p(\sigma^2|z^{(s+1)},y,\boldsymbol{X})\)を生成する
5. \(\boldsymbol{\beta}^{2(s+1)}\sim p(\boldsymbol{X}|z^{(s+1)},\sigma^{2(s+1)},y,\boldsymbol{X})\)を生成する

糖尿病の例

## Don't run it again if you've already run it
runmcmc<-!any(dir("source") %in% "diabetesBMA.RData")

if(!runmcmc){ load("source/diabetesBMA.RData") }

if(runmcmc){
  
  BETA<-Z<-matrix(NA,S,p)
  z<-rep(1,dim(X)[2] )
  lpy.c<-lpy.X(y,X[,z==1,drop=FALSE])
  for(s in 1:S)
  {
    for(j in sample(1:p))
    {
      zp<-z ; zp[j]<-1-zp[j]
      lpy.p<-lpy.X(y,X[,zp==1,drop=FALSE])
      r<- (lpy.p - lpy.c)*(-1)^(zp[j]==0)
      z[j]<-rbinom(1,1,1/(1+exp(-r)))
      if(z[j]==zp[j]) {lpy.c<-lpy.p}
    }
    
    beta<-z
    if(sum(z)>0){beta[z==1]<-lm.gprior(y,X[,z==1,drop=FALSE],S=1)$beta }
    Z[s,]<-z
    BETA[s,]<-beta
    if(s%%10==0)
    { 
      bpm<-apply(BETA[1:s,],2,mean) ; plot(bpm)
      cat(s,mean(z), mean( (y.te-X.te%*%bpm)^2),"\n")
      Zcp<- apply(Z[1:s,,drop=FALSE],2,cumsum)/(1:s)
      plot(c(1,s),range(Zcp),type="n") ; apply(Zcp,2,lines)
    }
  } 
  save(BETA,Z,file="source/diabetesBMA.RData")
}

説明変数の事後確率

\(\hat{\beta}_{bma}\)

• \(\beta\)の事後平均\(\hat{\beta}_{bma}=\sum_{s=1}^S\beta^{(s)}/S\)で近似
  • \(\beta\)の(ベイズ)モデル平均(Bayesian model averaged)推定量と呼ばれる
• 異なるzの値から得られる回帰パラメータの平均
  • 単一のモデルで推定した\(\beta\)よりも性能が良くなることが多い

beta.bma<-apply(BETA,2,mean,na.rm=TRUE)
print(beta.bma)

 [1]  0.0002349689 -0.0859065610  0.3423203546  0.1784059483 -0.3762771126
 [6]  0.2781232114 -0.0044069481  0.0266849669  0.4270552871  0.0056635027
[11]  0.0042001884  0.0146388944  0.0047480761 -0.0021133196 -0.0033337196
[16] -0.0003881057  0.0114732862  0.0132863289  0.0114692165  0.1015032151
[21]  0.0010538499  0.0094550788 -0.0014360247 -0.0241805632  0.0027968474
[26]  0.0025953718  0.0450453591  0.0082762134  0.0149977516  0.0133097496
[31]  0.0015971870 -0.0010502676  0.0025582091 -0.0001147017  0.0010068874
[36]  0.0061732992  0.0159645414 -0.0012871720  0.0003351759 -0.0012023777
[41]  0.0014035074  0.0011455487  0.0201005590  0.0051882624  0.0033907763
[46]  0.0005978504  0.0022609082  0.0111657754  0.0001771299 -0.0005611542
[51]  0.0051321670 -0.0210297438 -0.0059695252  0.0059841046 -0.0059704014
[56]  0.0171164352  0.0053458132  0.0054515426 -0.0031012637  0.0058710738
[61]  0.0008870141 -0.0180022426  0.0044849865  0.0040386379

予測値と実測値の平均二乗誤差

[1] 0.4825364

• \(Y,x\)が無関係な状況では?
• \(\boldsymbol{X}\)とは独立に順番を入れ替えて作ったベクトル\(\tilde{y}\)
• 変数減少法: 18個の説明変数が残った
• \(\mathrm{Pr}(z_j=1|\boldsymbol{y,X})=\sum z_l^{(s)}/S\)だと全ての説明変数で事後確率は1/2未満
• 変数減少法より、本来ない効果あると見積もることはない(らしい)